Effect of hydrophobic solutes on the hquid-hquid critical point 
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Jagla ramp particles, interacting through a ramp potential with two characteristic length scales, are 
known to show in their bulk phase thermodynamic and dynamic anomalies, similar to what is found 
in water. Jagla particles also exhibit a line of phase transitions separating a low density liquid phase 
and a high density liquid phase, terminating in a liquid-liquid critical point in a region of the phase 
diagram that can be studied by simulations. Employing molecular dynamics computer simulations, 
we study the thermodynamics and the dynamics of solutions of hard spheres (HS) in a solvent 
formed by Jagla ramp particles. We consider the cases of HS mole fraction xhs ~ 0.10,0.15 and 
0.20, and also the case xhs ~ 0.50 (a 1:1 mixture of HS and Jagla particles). We find a liquid-liquid 
critical point, up to the highest HS mole fraction; its position shifts to higher pressures and lower 
temperatures upon increasing xhs- We also find that the diffusion coefficient anomalies appear to 
be preserved for all the mole fractions studied. 

PACS numbers: 64.70. Ja,65.20.-w,66. 10. C- 



I. INTRODUCTION 



Liquid water exhibits highly unusual thermodynamic 
and dynamic behavior [l|-|3|. Among its most known 
anomalies there are the decrease in density upon isobaric 
cooling (density anomaly), the apparent divergences of 
thermodynamic response functions such as the isother- 
mal compressibility, the coefficient of thermal expansion 
and the isobaric specific heat upon cooling and the in- 
crease in diffusivity upon isothermal compression ( "diffu- 
sion anomaly" ) . It has been hypothesized that the ther- 
modynamic anomalies of water may be related to the 
presence of liquid-liquid (LL) critical point (LLCP) in 
the deeply supercooled region [i-^. This hypothesized 
LLCP is the end point of the liquid-liquid coexistence 
line separating two distinct liquid phases: a low density 
liquid (LDL) and a high density liquid (HDL). 

Molecular dynamics simulations of water aiming to 
study the LLCP and phenomena related to it often 
employ water models that reproduce the tetrahedral 
orientation-dependent interactions of real water. How- 
ever several papers have recently shown that tetrahedral- 
ity and even orientation-dependent interactions are not 
necessary conditions for the appearance of the density 
and the diffusion anomalies or the presence of a LL tran- 
sition [17T426]. There exists a family of spherically sym- 
metric potentials composed by a hard core and a linear 
repulsive ramp, the Jagla ramp model (27l [28| that can 
be tuned, varying the ratio between its two characteris- 
tic length scales, in order to span the range of behavior 
from hard spheres to water- like [l^, HI, 113 • With the ap- 
propriate choice of parameters the Jagla ramp potential 
with two characteristic length scales displays both ther- 
modynamic and dynamic anomalies and a LL transition. 

Spherically symmetric potentials with softened core 
have been used as coarse-grained models for a variety 
of substances beside water, such as metallic systems and 



colloidal suspensions [3lH33|. 

Besides its rather unique behavior as a pure liquid, wa- 
ter is also a remarkable solvent. In particular, phenomena 
related to the solvation of apolar solutes in water are in- 
teresting since they encompass biological membrane for- 
mation, globular protein folding and also the stability of 
mesoscopic assembly (s^ - ISTj . A large number of papers 
have in the past addressed the phenomenon of hydropho- 
bic hydration (see for example Refs. ItI. IbsI I38l - l49. ) . Small 
apolar solutes, such as alkanes or noble gases are poorly 
soluble in water. The solvation free energy of this kind 
of solutes is large and positive due to the large and nega- 
tive entropy contribution, the latter having been related 
to the structure of the hydrophobic hydration shell [i^ . 
These quantities have a marked temperature dependence 
and one intriguing anomaly of water in solutions is repre- 
sented by the increase in solubility of hydrophobic gases 
upon decreasing temperature [soj - 




FIG. 1. Spherically symmetric Jagla ramp potential. The 
potential has two scales: the hard core diameter r = a and 
the soft-core diameter r = h. In this case Ur/Uo ~ 3.56, 
b/a — 1.72 and c/a = 3. We have discretized the potential 
using a discretization step AU=Uo/8. 
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FIG. 2. (Color online) Isochores in the P — T plane for the four solutions at different hard spheres mole fractions, (a) 
xhs = 0.10, isochores are drawn for p = No? /V^ where L/a = VTA, 17.3, . . . , 15.5. The range of corresponding densities span 
from p — 0.328 to p = 0.464 (from bottom to top). The temperatures range is 0.255 < T < 0.380. (b) xhs = 0.15, isochores 
are drawn for p = No? /Lp' where L/a = 17.4, 17.3, . . . , 15.5. The range of corresponding densities span from p = 0.328 to 
p = 0.464 (from bottom to top). The temperatures range is 0.275 < T < 0.360. (c) xhs ~ 0.20, isochores are drawn for 
p = Na^/L^ where L/a = 17 A, 17.3, . . . , 15.5. The range of corresponding densities span from p — 0.328 to p = 0.464 (from 
bottom to top). The temperatures range is 0.275 <T< 0.360. (d) xhs = 0.50, isochores are drawn for densities p = Na"^ /L^ 
where L/a — 15.1, 15,0, ... , 13.7. The range of corresponding densities span from p — 0.502 to p = 0.672 (from bottom to 
top). The temperatures range is 0.210 < T < 0.300. In all panels the lines are fourth degree polynomial fits to simulated state 
points. For all mole fractions the position of the LLCP (circles), the LDL LMS (triangles up) and the HDL LMS (triangles 
down) are also reported. 



Experimentally the solubility of small apolar hy- 
drophobic solutes decreases upon decreasing temperature 
until a minimum is reached in the temperature range 
from 310 K to 350 K. Upon further teniperature decrease, 
the solubility increase monotonically [s^, • The solu- 
bility of model hydrophobic solutes, hard spheres, in the 
two-scale ramp potential particles has been recently as- 
sessed [S^ . It was found that the mixture of ramp parti- 
cles and hard spheres shows a temperature of minimum 
solubility, similarly to experimental results in water. The 
increased solubility upon cooling is connected to the pres- 
ence of two repulsive length scales in the model potential, 
the hard core corresponding to nearest-neighbor shell of 
solvent molecules and the soft repulsive core. Moreover 
it was observed that hard spheres are more favorably sol- 
vated in low density phases, in accord with what found 
in simulations of water 

The nature of critical phenomena in the presence of 
solutes has been extensively studied in literature with 



regard to the liquid-gas critical point [5ll - [53j |. while the 
effect of solutes on the LLCP of the solvent is a rela- 
tively new subject [s^. [sBj. In this work we investigate 
the thermodynamic and dynamic properties of the mix- 
ture of ramp potential supplemented by an attractive tail 
and apolar solutes modeled by hard spheres. We exam- 
ine three mole fractions of hard spheres, xhs = 0.10, 0.15 
and 0.20. We also investigate thermodynamics and dif- 
fusivity for a 1:1 mixture of hard spheres and Jagla ramp 
particles, xhs = 0.50. The paper is structured as follows. 
In Sec. HI] we give the details of the interaction potential 
and of the simulation method. We present results and 
discussion in Sec. IIIII and conclusions in Sec. |IV] 



II. METHODS 

We perform discrete molecular dynamics [2l|, [H, 
on systems composed by iV = Nramp + Nhs = 1728 
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FIG. 3. (Color online) Isotherms in the P ~ p plane for the four solutions at different hard spheres mole fractions, (a) 
XHS = 0.10, isotherms are drawn from T = 0.255 to T = 0.380 every AT = 0.005 (from bottom to top) . (b) xhs = 0.15, 
isotherms are drawn from T — 0.275 to T = 0.360 every AT = 0.005 (from bottom to top), (c) xhs = 0.20, isotherms are 
drawn from T = 0.275 to T = 0.360 every AT — 0.005 (from bottom to top), (d) xhs ~ 0.50, isotherms are drawn from 
T = 0.210 to T = 0.300 every AT = 0.005 (from bottom to top). In all panels the lines are fourth degree polynomial fits 
to simulated state points. The circles represent the position of the LLCP. Every other line has been made dashed to help 
distinguishing in between them. 



particle, where Nramp is the number of water-like parti- 
cles and Nhs is the number of hard spheres. We study 
four systems with different composition. The solute con- 
tent in the four systems is Nhs = 173 for x^s = 0.10, 
Nhs = 260 for xhs = 0.15, Nhs = 345 for xhs = 0.20 
and Nhs = 864 for xhs = 0.50 

The pair- wise Jagla ramp interaction potential [13, [1^ 
has two characteristic length scales, the hard-core dis- 
tance r = a and the soft-core distance r = b. The min- 
imum of the energy Uo is corresponds to the soft-core 
distance. An attractive tail extends up to r = c. The 
potential has been discretized, in order to be able to em- 
ploy the algorithm of discrete molecular dynamics. We 
partition the repulsive ramp into 36 steps of width 0.02a 
and the attractive ramp into eight steps of width 0.16a. 
The AU at each step is Uo/8 = 0.125. The parameters of 
the ramp potential |21|] have been set to b/a = 1.72 and 
c/a = 3. Ub. = 3.56 J7o is defined as the value of the en- 
ergy at r = a, obtained via least-squares linear fit to the 
discretized repulsive ramp. The shape of the spherically 
symmetric Jagla ramp potential employed in this work 
is shown in Fig. [TJ As in previous papers [13, [HI , this 



parametrization of the ramp potential prevents the oc- 
currence of crystallization. Furthermore, with this choice 
of parameters, the line of LL phase coexistence extends 
into the equilibrium liquid phase and ends in a LLCP. 
The diameter of the hard spheres is a, the same as the 
hard-core distance of Jagla ramp interaction potential. 
The solvent and the solute interact via a hard-core po- 
tential. 

We express all quantities in reduced units. Distances 
are in unit of a, energies in units of Uq and time in units 
a-^/m/C/o, where m is the mass that is assumed to be 
unitary. The density defined as p = N/L^, where L is 
the edge of the cubic simulation box, is measured in units 
of a~^, pressure in unit of Uq/o?' and temperature in units 
of Uo/kB- 

We perform simulations at constant N, V and T, with 
T controlled by rescaling the velocity of the particles 
with a modified Berendsen algorithm (details are given 
in Ref . f56h . or at constant N, P and T, with P controlled 
by allowing the edge of the simulation box to vary with 
time applying the standard Berendsen algorithm [56[. 
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FIG. 4. Position in the P — T plane of tiie LLCP (squares) 
and of LDL (triangles up) and HDL (triangles down) LMS 
lines for the solutions of hard spheres in Jagla ramp particles, 
with XHS = 0.10,0.15 and 0.20. The position of the LLCP 
of bulk Jagla ramp particles is also reported for comparison. 
In the inset the analogous quantities are shown for the 1:1 
mixture of Jagla ramp particles and hard spheres. 



III. RESULTS AND DISCUSSION 

We study the thermodynamics of the solutions of hard 
spheres in Jagla ramp potential water-like particles, ana- 
lyzing the isochores in the P — T plane and the isotherms 
in the P — p plane. We study ranges of temperature and 
density for which the hard spheres are completely soluble 
in the Jagla ramp potential liquid [3^ . As a consequence 
in the range we study here, no solvent-solute demixing 
occurs. The position of the LLCP has been estimated 
considering the inflection point in the isotherms where 



dP 
'dp 



92p 



= 



(1) 



The points of the LL limit of mechanical stability (LMS) 
lines have been determined by the points for which 
(dP/dp)T = 0, and {d'^P/dp'^)T 7^ where the isother- 
mal compressibility Kt = {dp/ dP)T / p diverges. 

The thermodynamic properties of the bulk Jagla ramp 
potential particles, with the same set of parameters used 
in this work, have been previously assessed (2^. In par- 
ticular the coordinates in the thermodynamic plane of 
the LLCP of bulk Jagla ramp potential particles are 
Tc = 0.375, Pc = 0.243 and pc = 0.37. 

In Fig. [5] we present the isochores P — T plane for all 
the solutions at different mole fractions. In this figure 
the location of the LLCP and the two branches of the 
LL LMS lines are also shown. For mole fractions up to 
xhs — 0.20 we observe that the isochores converge very 
clearly to the LLCP and cross at points corresponding to 
the LMS lines. For mole fractions xhs = 0.10,0.15 and 
0.20 we also note that the low density isochores appear 



TABLE L Position of the LLCP for bulk Jagla ramp potential 
particles and for the solutions of hard spheres in Jagla ramp 
potential particles with xhs = 0.10, 0.15 and 0.20 and for 1:1 
mixture of Jagla ramp potential particles and hard spheres 
{xhs = 0.50). 



Xhs 




Pc 


Pc 


0.00 (bulk) 


0.375 


0.243 


0.370 


0.10 


0.346 
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0.320 


0.362 
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0.230 


0.645 


0.597 



almost flat, signaling the presence of the density anomaly. 
The points of the temperature of maximum density line 
are in fact given by the zeros of the coefficient of thermal 
expansion, where {dP/dT)p — 0. In the case of the mix- 
ture with Xhs = 0.50 the convergence of the isochores to 
the LLCP appears less precise. In this case we observe 
that the isochores are not flat any longer, indicating the 
disappearance of the density anomaly at this high mole 
fraction of hard spheres. 

Fig. [3] shows the isotherms of the mixtures with xhs = 
0.10,0.15,0.20 and 0.50 along with the position of the 
LLCP in the P ~ p plane. We can observe the inflec- 
tion point in the isotherms corresponding to the critical 
point below which the isotherms show van der Waals- 
like loops indicating LL coexistence. We also observe 
that upon increasing the mole fraction of solutes the den- 
sity range of the region of coexistence narrows. We also 
point out that for low density the isotherms cross due 
to density anomaly for mole fractions xhs = 0.10,0.15 
and 0.20 . In fact for the points where isotherms cross, 
{dP/dT)p = 0, so the crossing of the isotherms indicates 
a density anomaly. Also in the case of xhs = 0.50 mix- 
ture the presence of a LLCP can be highlighted by the 
inflection point in the critical isotherm and the van der 
Waals-like loop structure of the isotherms below the crit- 
ical temperature. An important feature of isotherms of 
this system is that they do not cross, confirming what we 
have found looking at the isochores in Fig. For this 
very high mole fraction of solutes the density anomaly is 
not longer present, while the LLCP is still found. 

The coordinates of the estimated positions of the 
LLCP for bulk Jagla ramp potential particles [l^l and for 
all the systems with different composition are reported 
together in Table I. 

The positions of the LLCP and the two branches of the 
LL LMS line for the solutions with xhs = 0.10, 0.15 and 
0.20 and for the 1:1 system are also reported in Fig. |4] 
in the P — T plane. The position of the LLCP moves to 
higher pressures and lower temperatures upon increas- 
ing the solute mole fraction. This shift of the critical 
point could be connected to the fact that hard spheres 
are more favorably solvated in LDL 0, [s^l . In fact the 
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FIG. 5. (Color online) Diffusion coefficients of Jagla ramp particles at constant temperature for the four solutions at different 
hard spheres mole fractions, (a) xhs ~ 0.10, isotherms of the diffusion coefficient are drawn from T = 0.255 to T = 0.380 every 
AT — 0.005 (from bottom to top), (b) xhs = 0.15, isotherms of the diffusion coefficient are drawn from T = 0.275 to T = 0.360 
every AT — 0.005 (from bottom to top), (c) xhs ~ 0.20, isotherms of the diffusion coefficient are drawn from T = 0.275 to 
T — 0.360 every AT = 0.005 (from bottom to top), (d) xhs = 0.50, isotherms of the diffusion coefficient are drawn from 
T = 0.210 to T = 0.300 every AT = 0.005 (from bottom to top). For all panels, the lines are fourth degree polynomial fits to 
simulated state points. The dashed lines join the diffusivity extrema (maxima for xhs = 0.10,0.15,0.50, maxima and minima 
for xhs = 0.20). 



LDL region of the phase diagram progressively widens 
upon increasing the mole fraction of hard spheres. 

Looking at the isochores (Fig. [2]) and at the isotherms 
(Fig. [3]) we can observe that upon increasing the mole 
fraction of hard spheres, the width of the coexistence 
envelope is reduced. In fact the region of crossing of 
the isochores is progressively reduced upon increasing the 
solute mole fraction, on going from the system at xhs = 
0.10 to the one at xhs = 0.50. Also the width of the 
loop region in the isotherms plane becomes more narrow 
spanning a minor range of densities, upon increasing the 
solute mole fraction. From these considerations we can 
argue that upon a further increase in solute mole fraction, 
the LL critical phenomenon will disappear. This is in 
agreement with studies of the LL transition in aqueous 
solutions (55l.l57l-[59t. 

In Fig. [S] we present the diffusion coefficients for the 
Jagla ramp potential particles, calculated at constant 
temperature as a function of the density of the system, 
for all the mole fractions studied. For all solute mole frac- 
tions, we find the appearance of the diffusivity anomaly. 



The Xhs = 0.10,0.15 and 0.20 mixtures are studied in 
the same set of densities, p = 0.328 to p = 0.464. In this 
range maxima are evident for xhs = 0.10,0.15 and 0.20, 
while minima can be seen only for the xhs = 0-20 mix- 
ture. Minima for the xhs = 0-10 and 0.15 mixtures are 
to be found at lower densities, out of the spanned density 
range. This is due to the narrowing of the region of the 
diffusivity anomaly upon increasing the solute mole frac- 
tion. Thus not only thermodynamic anomalies but also 
dynamic anomalies exist in a smaller range of densities, 
upon increasing the solute mole fraction. The diffusion 
coefficients at constant pressure for the xhs — 0.50 are 
studied in the density range p = 0.502 to p = 0.672, thus 
cannot be directly compared with mixtures with lower 
solute mole fraction as the range of densities spanned 
is different. However we can see that they also exhibit 
diffusion anomaly with the presence of maxima in the 
isotherms of the diffusion coefficient. 

In Fig. ini we report the behavior of the diffusion co- 
efficient of Jagla ramp particles and hard spheres cal- 
culated for a constant pressure path above the critical 
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FIG. 6. Diffusion coefScients for hard spheres and Jagla ramp 
particles along a constant pressure path above the critical 
point. P = 0.321 for xhs = 0.10, P = 0.347 for xhs = 0.15 
and P = 0.382 for xhs ~ 0.20. Lines correspond to the 
Arrhenius fit for the Jagla ramp particles. In this figure the 
diffusion coefficients at constant pressure for the bulk Jagla 
ramp particles at P = 0.250 and P = 0.275 are also reported 
for comparison. 



pressure for the solutions with hard sphere mole fractions 
Xhs = 0.10, 0.15 and 0.20. For all compositions the pres- 
sure is set to the critical pressure plus AP — 0.020. In 
the bulk Jagla ramp particles system it was found that 
the trend of the diffusion coefficient, calculated on a cool- 
ing path above the critical point, shows a crossover from 
a high temperature behavior (LDL-like) to a stronger 
(HDL-like) behavior [2^. This change in trend was con- 
nected to the maximum of the specific heat that oc- 
curs at the Widom line via the Adam-Gibbs equation 
D = Do exp(— C/T5con/), where Sconf is the config- 
urational entropy. We can observe that the crossover 
in the behavior of the diffusion coefficient is maintained 
up to the highest mole fraction studied {xhs = 0.20), 
and at low enough temperature, below the Widom line, 
for all compositions the diffusion behavior becomes that 
of a HDL Jagla liquid. The temperature at which the 
dynamic crossover occurs decreases upon increasing the 
mole fraction of solutes. As the Widom line is connected 
to the LLCP, the shift to lower temperatures of the dy- 
namic crossover confirms the shift to lower temperatures 
of the LLCP we have found studying the thermodynam- 
ics (see Fig. U). The trend of the diffusion coefficient 
of hard spheres closely follows that of the Jagla ramp 
potential particles, thus we can derive that the diffusive 
behavior of the hard spheres solute is determined by the 
solvent. 

Therefore we conclude that the dynamic cross-over 
found in bulk Jagla ramp particles system (analogous 
to liquid water) at the Widom line, is preserved in the 
solutions and that the temperature of dynamic cross-over 
decreases, upon increasing the mole fraction of hydropho- 
bic solutes. 



Finally the solute-solute radial distribution functions 
are shown in Fig. [T] The gHS-Hsii^) have been re- 
ported for the critical density of the systems and for 
three temperatures, T=0.360, the critical temperature 
and T=0.290 for the solutions with xhs 0.10,0.15 
and 0.20 and for T=0.300, the critical temperature and 
T=0.210 for Xhs — 0.50. A small tendency of the so- 
lutes to cluster upon increasing concentration can be seen 
looking at the progressive increase in the g{r) near the 
hard-core distance. At the lowest temperature the solute- 
solute radial distribution functions also show anomalous 
behavior with a progressive decay toward the asymptotic 
value of 1 at large r that tends to disappear upon increas- 
ing the mole fraction of solutes. The decay indicates seg- 
regation of the HS into the LDL phase below the critical 
point. 



IV. CONCLUSIONS 

We performed discrete molecular dynamics simulations 
on the mixture of Jagla ramp potential (water-like) par- 
ticles and hard spheres, at four different compositions, 
XHS = 0.10,0.15,0.20 and 0.50. The thermodynamic 
and dynamic behavior were studied for the solutions with 
Xhs — 0.10,0.15 and 0.20. Thermodynamics and diffu- 
sive behavior were also studied for the 1:1 mixture of 
Jagla ramp particles and hard spheres. 

The analysis of the isochores and the isotherms plane 
revealed the presence of a liquid-liquid critical point 
for all the investigated system. We also found that 
while density anomaly is present in the solutions with 
Xhs = 0.10,0.15 and 0.20, it disappears for the 1:1 mix- 
ture. Furthermore we also observe a narrowing of the 
coexistence envelope in both planes, upon increasing the 
solute mole fraction. The position of the critical point 
is found to shift toward lower temperatures and higher 
pressures, upon increasing the hard sphere mole fraction. 
This shift may be related to the favored solvation of hard 
spheres in the LDL Jagla ramp potential particles sol- 
vent. 

The complete phase diagram of the solution in the 
vicinity of the LLCP may be quite complex and might 
include also an upper critical solution temperature dif- 
ferent from the LLCP. The interplay of these two critical 
points is an interesting subject which deserves separate 
investigation. 

The appearance of extrema in the behavior of the con- 
stant temperature diffusion coefficient for Jagla ramp po- 
tential particles, i.e. the diffusion anomaly is also pre- 
served up to the highest mole fraction of hard spheres 
studied, with a narrowing of the anomalous region upon 
increasing the solute mole fraction, also for this dynamic 
property. 

Finally a change in trend of the constant pressure diffu- 
sion coefficient from LDL-like behavior to HDL-like be- 
havior can be observed when cooling the system in a 
constant pressure path, above the critical pressure. This 
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FIG. 7. Hard sphere-hard sphere radial distribution functions for the mixtures with xhs = 0.10,0.15,0.20 and 0.50 calculated 
at their critical densities (see Table I). The QHS-Hsif) are reported at temperatures T=0.360, 0.345 and 0.290 for xhs ~ 0.10, 
T=0.360, 0.330 and 0.290 for xhs = 0.15, T=0.360, 0.320 and 0.290 for xhs = 0.15 and T=0.300, 0.230 and 0.210 for 
XHS = 0.50. 



cross-over in the dynamic behavior can be related to the 
crossing of the Widom line, above the critical point. The 
dynamic cross-over observed for bulk Jagla ramp parti- 
cles at the Widom line shifts to lower temperatures upon 
increasing the content of solutes. 
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